Epidemiology and Clinical Outcomes of Patients with Acute Hypoxaemic Respiratory Failure Admitted to Australian and New Zealand Intensive Care Units
Statistical Analysis
Author
Dr Benjamin Moran, MBBS, BMedSci (Hons), MMedStats, FCICM, FANZCA
Published
June 21, 2024
1 Introduction
This is an explanation of the statistical analysis for the study examining the epidemiology and outcomes of patients admitted to ICU with acute hypoxaemic respiratory failure (AHRF).
2 Methods
This is a retrospective study using data from the Australian and New Zealand Intensive Care Society (ANZICS) adult ICU patient database (APD). This manuscript has been prepared and reported in accordance with the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) Statement.
2.1 Patient Population
We included all patients in the ANZICS APD from 01/01/2005 to 31/12/2022. Patients were categorised into PF ratio subgroups according to the severity definitions below.
2.2 Aims and Objectives
Determine the association between PaO2:FiO2 ratio and hospital mortality.
Determine the hospital mortality trends over time across PaO2:FiO2 ratio categories.
Determine the admissions to ICU with AHRF over time across PaO2:FiO2 ratio categories.
Determine the trends in proportion of hospital outcomes (Hospital mortality, ICU mortality, ICU length of stay, hospital length of stay, discharge home, nursing home/palliation, rehabilitation, other acute care hospital, other hospital ICU) over time by AHRF categories.
Determine the association between PaO2:FiO2 ratio and hospital mortality in pre-specified subgroups (Level of ventilatory support, gender, age categories, country, hospital type, admission diagnosis, presence of treatment limitations at ICU admission)
Determine the time to death of each AHRF categories.
Determine the time to death of time categories.
2.3 Definition of Acute Hypoxaemic Respiratory Failure
Currently, there is heterogeneity in the precise definition of Acute Hypoxaemic Respiratory Failure (AHRF). For this study, we have defined AHRF as a PaO2:FiO2 ratio < 300 mmHg. AHRF can further be categorised as mild (PaO2:FiO2 200-300 mmHg), moderate (PaO2:FiO2 100-200 mmHg) and severe (PaO2:FiO2 <100 mmHg). Further categorsiation can be made into moderate-severe (PaO2:FiO2 100-150 mmHg).
2.4 Statistical Analysis
Baseline ICU- and patient-level characteristics and unadjusted outcomes were summarized using standard descriptive statistics. Continuous variables were reported as either means with standard deviation or medians and interquartile range, and categorical variables reported as number and percentages.
2.4.1 Covariate Selection for Multivariable Regression Models for Adjusted Outcomes
A directed acyclic graph (DAG) was used to generate a minimum adjustment set of covariates for the causal pathway from AHRF to hospital outcomes. Covariates were selected based on potential mechanistic associations with other variables. In addition to this, variables were selected if they increased the precision of the estimate. All back-door paths were closed and a minumum adjustment set of selected variables were used for that particular model. For all outcomes, baseline variables that were included in the model included PaO2:FiO2 ratio, chronic respiratory disease, chronic cardiovascular disease, frailty, smoking intensity and severity of illness scores.
2.4.2 Model Fitting
To analyse the association of AHRF and hospital outcomes, a multivariable, hierarchical logistic regression model was used with patients nested within sites and site treated as a random effect. Predicted probabilities were generated from the model output and displayed graphically as the effect on hospital outcomes by either the continuous PaO2:FiO2 ratio or time. The association between either PaO2:FiO2 ratio or time and hospital outcomes was modeled using restricted cubic splines with 4 knots to allow for non-linear association. Changes over time were described with the estimand of Absolute Risk Reduction (with 95% CI). To determine whether changes over time differed between different AHRF categories, an interaction term between time and the AHRF categories was fitted.
2.4.3 Time-to-Event Analysis
Time to death was analysed using Cox-proportional hazards regression with random effects (frailty model), with covariates included as determined by the minimum adjustment set from the previously constructed DAG. The results are presented as a Kaplan-Meier curve generated from the Cox Proportional Hazards regression with the HR (95% CI) for each AHRF category compared to no AHRF (“None” category). To further investigate the change over time, the study period was divided into 3 cohorts (2005-09, 2010-14, 2015-19, 2020-2022). These were presented as Kaplan-Meier curves from the Cox regression, with the HR (95% CI) comparing each time epoch to the initial time (2005-09). Each AHRF category is displayed as a separate Kaplan-Meier curve to examine the change within that category over time.
As there were >1,500,000 patients in the dataset, a 2-sided p-value of 0.001 was used for statistical significance. Given that there is an increased risk of Type-1 error with multiple testing, the results of the secondary objectives should be viewed as exploratory. Hence, no adjustment for multiplicity was used. Only patients with complete data for all covariates were included in the analysis. Statistical analyses were performed using R Version 4.3.1 (R Core Team, R Foundation for Statistical Computing, Vienna, Austria) and RStudio Version 2023.06.1 (Posit Software, PBC, Boston, MA). Packages used for analysis included tidyverse, data.table, gtsummary, lme4, survival and ggsurvfit.
2.5 Subgroup Analysis
Patients were analysed for the association of PF ratio and in-hospital mortality in the following subgroups: Receiving invasive ventilation on day 1, receiving invasive ventilation at any time during the ICU admission, receiving ECMO during the ICU admission, levels of ventilatory support (none, non-invasive ventilation, invasive ventilation, extra-corporeal membrane oxygenation), gender, age categories, admission diagnosis (medical, surgical, cardiac, neurosurgery/trauma, COVID pneumonitis), and treatment limitation status on ICU admission.
2.6 Sensitivity Analysis
Post-hoc sensitivity analyses were performed on the pre-specified variables of severity of illness scores, and for duration of invasive ventilation. In the first sensitivity analysis, modelling was repeated substituting the APACHE III score for APACHE II, ANZROD and SOFA scores. In the second sensitivity analysis, modelling was repeated on the pre-specified variables of invasive ventilation of more than 12 hours.
3 Results
Of 1,560,221 patients admitted to 211 ICUs during the study period, 826,106 (52.9%) patients had acute hypoaxemic respiratory failure. Of this cohort, 424382 (27.2%) had mild AHRF, 324,630 (20.8%) had moderate AHRF and 77,094 (4.9%) had severe AHRF.
3.1 Patient Demographics
Below are the demographic tables. This table has the 4 AHRF categories (none, mild, moderate, severe) to look at the breakdown of patients within each AHRF category.
arf |>group_by(IcuAdmitYYYY, ahrf_cause) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(!ahrf_cause=="NA") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = ahrf_cause))+geom_line()+labs(x="Year", y="Proportion of ICU Admission (%)", title="ICU Admission Diagnoses over Time", color ="ICU Admission Diagnosis")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 50, by =5))+ylim(0,50)
`summarise()` has grouped output by 'IcuAdmitYYYY'. You can override using the
`.groups` argument.
Scale for y is already present. Adding another scale for y, which will replace
the existing scale.
arf |>group_by(IcuAdmitYYYY, pf_cat, ahrf_cause) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(ahrf_cause=="Cardiac Surgery") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Proportion of Disease-Specific Admission (%)", title="ICU Admission After Cardiac Surgery over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 60, by =5))+ylim(0,60)
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
Scale for y is already present. Adding another scale for y, which will replace
the existing scale.
arf |>group_by(IcuAdmitYYYY, pf_cat, ahrf_cause) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(ahrf_cause=="Respiratory Disease") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Proportion of Disease-Specific Admission (%)", title="ICU Admission With Respiratory Disease over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 35, by =5))+ylim(0,35)
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
Scale for y is already present. Adding another scale for y, which will replace
the existing scale.
arf |>group_by(IcuAdmitYYYY, pf_cat, ahrf_cause) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(ahrf_cause=="Cardiovascular Disease") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Proportion of Disease-Specific Admission (%)", title="ICU Admission With Cardiovascular Disease over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 30, by =5))+ylim(0,30)
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
Scale for y is already present. Adding another scale for y, which will replace
the existing scale.
arf |>group_by(IcuAdmitYYYY, pf_cat, ahrf_cause) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(ahrf_cause=="Sepsis") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Proportion of Disease-Specific Admission (%)", title="ICU Admission With Sepsis over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 20, by =5))+ylim(0,20)
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
Scale for y is already present. Adding another scale for y, which will replace
the existing scale.
arf |>group_by(IcuAdmitYYYY, pf_cat, ahrf_cause) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(ahrf_cause=="Gastrointestinal Surgery") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Proportion of Disease-Specific Admission (%)", title="ICU Admission After Gastrointestinal Surgery over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 20, by =5))+ylim(0,20)
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
Scale for y is already present. Adding another scale for y, which will replace
the existing scale.
arf |>group_by(IcuAdmitYYYY, pf_cat, ahrf_cause) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(ahrf_cause=="Other Non-Operative") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Proportion of Disease-Specific Admission (%)", title="ICU Admission With Other Non-Operative over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 30, by =5))+ylim(0,30)
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
Scale for y is already present. Adding another scale for y, which will replace
the existing scale.
arf |>group_by(IcuAdmitYYYY, pf_cat, ahrf_cause) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(ahrf_cause=="Other Post-Operative") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Proportion of Disease-Specific Admission (%)", title="ICU Admission With Other Post-Operative over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 40, by =5))+ylim(0,40)
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
Scale for y is already present. Adding another scale for y, which will replace
the existing scale.
3.4 Admissions to ICU with AHRF Over Time
The proportion of AHRF admissions to ICU was 53.6% (559/1,042) in 2005 and 52.0% (66521/127,934) in 2022. The proportion of patients admitted with mild AHRF to ICU was 27.2% (283/1,042) in 2005 and 27.6% (35,306/127,934) in 2022. The proportion of patients admitted with moderate AHRF to ICU was 20.6% (215/1,042) in 2005 and 19.9% (25,508/127,934) in 2022. The proportion of patients admitted with severe AHRF to ICU was 5.9% (61/1,042) in 2005 and 4.5% (5,707/127,934) in 2022.
arf |>mutate(ahrf =if_else(pf_cat =="None", 0,1)) |>group_by(IcuAdmitYYYY,ahrf) |>summarise(n_ahrf =n()) |>mutate(percent =100*n_ahrf/sum(n_ahrf)) |>filter(ahrf=="1") |>ggplot()+geom_line(aes(x = IcuAdmitYYYY, y = percent))+labs(x="Year", y="ICU Admission (%)", title="Proportion of Patients Admitted to ICU with AHRF over Time")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+coord_cartesian(ylim =c(0, 100))+scale_y_continuous(breaks =seq(0, 100, by =10))
`summarise()` has grouped output by 'IcuAdmitYYYY'. You can override using the
`.groups` argument.
arf |>group_by(IcuAdmitYYYY,pf_cat) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="ICU Admission (%)", title="Proportion of Patients Admitted to ICU over Time", subtitle ="By AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+coord_cartesian(ylim =c(0, 100))+scale_y_continuous(breaks =seq(0, 100, by =10))
`summarise()` has grouped output by 'IcuAdmitYYYY'. You can override using the
`.groups` argument.
3.5 Unadjusted Outcomes
3.5.1 All Outcomes
Below are the unadjusted outcomes for the different AHRF categories. The first table contains the outcomes for patients with AHRF compared to those without it, and the second table has the 4 categories (none, mild, moderate, severe).
3.5.2 Unadjusted In-Hospital Mortality and Continuous PaO2:FiO2 Ratio
The figures below represent the unadjusted outcomes.
ggplot(data=arf, aes(x=pf_cont, y=DIED_HOSP))+geom_smooth(data = arf, method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Hospital Mortality (%)", title="Unadjusted Association of Day 1 PF Ratio and Hospital Mortality")+coord_cartesian(ylim =c(0,0.45))+scale_y_continuous(labels = scales::percent)+theme_bw() +theme(panel.grid.major.x =element_blank())
3.5.3 Unadjusted Trends in Hospital Outcomes Over Time by AHRF Categories
ggplot(data=arf, aes(x=IcuAdmitYYYY, y=DIED_HOSP))+geom_smooth(data = arf, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf, aes(x=IcuAdmitYYYY, y=DIED_HOSP, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Hospital Mortality (%)", title="Unadjusted Hospital Mortality over Time", color ="AHRF Category")+scale_y_continuous(labels = scales::percent)+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))
ggplot(data=arf, aes(x=IcuAdmitYYYY, y=DIED_ICU))+geom_smooth(data = arf, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf, aes(x=IcuAdmitYYYY, y=DIED_ICU, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="ICU Mortality (%)", title="Unadjusted ICU Mortality over Time", color ="AHRF Category")+scale_y_continuous(labels = scales::percent)+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))
ggplot(data=arf, aes(x=IcuAdmitYYYY, y=icu_los))+geom_smooth(data = arf, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf, aes(x=IcuAdmitYYYY, y=icu_los, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="ICU LOS (Days)", title="Unadjusted ICU LOS over Time", color ="AHRF Category")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))
ggplot(data=arf, aes(x=IcuAdmitYYYY, y=hosp_los))+geom_smooth(data = arf, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf, aes(x=IcuAdmitYYYY, y=hosp_los, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Hospital LOS (Days)", title="Unadjusted Hospital LOS over Time", color ="AHRF Category")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))
arf |>group_by(IcuAdmitYYYY, hosp_dc_dest) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(!hosp_dc_dest=="Died") |>filter(!hosp_dc_dest=="Mental Health") |>filter(!hosp_dc_dest=="Hospital in the Home") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = hosp_dc_dest))+geom_line()+labs(x="Year", y="Hospital Outcome (%)", title="Hospital Discharge Destination over Time", color ="Hospital Discharge Destination")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 70, by =5))
`summarise()` has grouped output by 'IcuAdmitYYYY'. You can override using the
`.groups` argument.
arf |>group_by(IcuAdmitYYYY, pf_cat, hosp_dc_dest) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(hosp_dc_dest=="Home") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Hospital Outcome (%)", title="Hospital Discharge to Home over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 80, by =5))
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
arf |>group_by(IcuAdmitYYYY, pf_cat, hosp_dc_dest) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(hosp_dc_dest=="Nursing Home/Chronic Care/Palliative Care/Rehabilitation") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Hospital Outcome (%)", title="Hospital Discharge to Nursing Home/Palliative/Rehabilitation Care over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 16, by =2))
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
arf |>group_by(IcuAdmitYYYY, pf_cat, hosp_dc_dest) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(hosp_dc_dest=="Other Acute Care Hospital") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Hospital Outcome (%)", title="Hospital Discharge to Other Hospital over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 30, by =5))
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
arf |>group_by(IcuAdmitYYYY, pf_cat, hosp_dc_dest) |>summarise(n =n()) |>mutate(percent =100*n/sum(n)) |>filter(hosp_dc_dest=="Other Hospital- ICU") |>ggplot(aes(x = IcuAdmitYYYY, y = percent, color = pf_cat))+geom_line()+labs(x="Year", y="Hospital Outcome (%)", title="Hospital Discharge to Other ICU over Time", subtitle ="by AHRF Category", color ="AHRF Category")+theme_bw() +theme(panel.grid.minor.x =element_blank(), panel.grid.minor.y =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 5.5, by =1))
`summarise()` has grouped output by 'IcuAdmitYYYY', 'pf_cat'. You can override
using the `.groups` argument.
3.6 Adjusted Outcomes
3.6.1 Primary Outcome: Conditional Probability of In-Hospital Mortality and AHRF Status
Adjusted mortality
3.6.2 Primary Outcome: Conditional Probability of In-Hospital Mortality and Continuous PaO2:FiO2 Ratio
Adjusted mortality increased as PaO2:FiO2 Ratio decreased (Adjusted ARR -3.598 (95%CI: -3.719 - -3.478)). This means that for every decrease of 100 in the PaO2:FiO2 Ratio, the conditional hospital mortality will increase by 3.598%. However, this effect is non-linear, as seen from the figure below.
Overall hospital mortality is decreasing over time. Over the study period, there was a reduction in adjusted hospital mortality from 13.278 (95% CI: 13.092 - 13.464) to 8.2 (95% CI: 8.128 - 8.271). This was an absolute risk reduction (ARR) of 5.079 (95% CI: 5.077 - 5.081). This change in outcome for each AHRF category over time were all significant (all p<0.001 for interactions compared to None category of AHRF), and the comparisons between each category and no-AHRF were all significant (p<0.001).
hosp_mort_time
# Remove NAsarf_icu_mort <- arf[complete.cases(arf$DIED_ICU),]# Run MLMarf_model_icu_mort <-glmer(as.numeric(DIED_ICU) ~0+ IcuAdmitYYYY + pf_cont + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),family =binomial(link ="logit"), data = arf_icu_mort, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))arf_model_icu_mort_tidy <-broom.mixed::tidy(arf_model_icu_mort, conf.int = T, exponentiate=TRUE)# Run MLM for pg_cat interaction termarf_model_icu_mort_cat <-glmer(as.numeric(DIED_ICU) ~0+ IcuAdmitYYYY:pf_cat + IcuAdmitYYYY + pf_cat + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),family =binomial(link ="logit"), data = arf_icu_mort, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))arf_model_icu_mort_cat_tidy <-broom.mixed::tidy(arf_model_icu_mort_cat, conf.int = T, exponentiate=TRUE)# Create predicted probabilities of hospital mortalityarf_icu_mort$pred_icu_mort <-predict(arf_model_icu_mort, type ="response")icu_mort_time <-ggplot(data=arf_icu_mort, aes(x=IcuAdmitYYYY, y=pred_icu_mort))+geom_smooth(data = arf_icu_mort, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf_icu_mort, aes(x=IcuAdmitYYYY, y=pred_icu_mort, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional Probability of ICU Mortality", title="ICU Mortality over Time", color ="AHRF Category")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))icu_mort_arr <-arr(icu_mort_time)
Overall ICU mortality is decreasing over time. Over the study period, there was a reduction in adjusted ICU mortality from 8.629 (95% CI: 8.472 - 8.786) to 5.538 (95% CI: 5.478 - 5.598). This was an absolute risk reduction (ARR) of 3.091 (95% CI: 3.089 - 3.093). This change in outcome between the categories of AHRF over time were all significant (all p<0.001 for interactions compared to None category of AHRF).
icu_mort_time
# Remove NAsarf_icu_los <- arf[complete.cases(arf$icu_los),]# Run MLMarf_model_icu_los <-lmer(icu_los ~0+ IcuAdmitYYYY + pf_cont + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),data = arf_icu_los)arf_model_icu_los_tidy <-broom::tidy(arf_model_icu_los, conf.int = T)# Run MLM for pg_cat interaction termarf_model_icu_los_cat <-lmer(icu_los ~0+ IcuAdmitYYYY:pf_cat + IcuAdmitYYYY + pf_cat + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),data = arf_icu_los)arf_model_icu_los_cat_tidy <-broom::tidy(arf_model_icu_los_cat, conf.int = T)# Create predicted probabilities of hospital mortalityarf_icu_los$pred_icu_los <-predict(arf_model_icu_los, type ="response")icu_los_time <-ggplot(data=arf_icu_los, aes(x=IcuAdmitYYYY, y=pred_icu_los))+geom_smooth(data = arf_icu_los, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf_icu_los, aes(x=IcuAdmitYYYY, y=pred_icu_los, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional ICU LOS (Days)", title="ICU Length of Stay over Time", color ="AHRF Category")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))icu_los_arr <-arr(icu_los_time)/100
Overall ICU length of stay is decreasing over time. Over the study period, there was a reduction in adjusted ICU LOS from 3.91562 (95% CI: 3.89709 - 3.93416) to 3.26592 (95% CI: 3.25883 - 3.27302). This was an absolute risk reduction (ARR) of 0.6497 (95% CI: 0.6495 - 0.6499). This change in outcome between the categories of AHRF over time were not significant (95% CI crossed 0 for interactions compared to None category of AHRF).
icu_los_time
# Remove NAsarf_hosp_los <- arf[complete.cases(arf$hosp_los),]# Run MLMarf_model_hosp_los <-lmer(hosp_los ~0+ IcuAdmitYYYY + pf_cont + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),data = arf_hosp_los)arf_model_hosp_los_tidy <-broom::tidy(arf_model_hosp_los, conf.int = T)# Run MLM for pg_cat interaction termarf_model_hosp_los_cat <-lmer(hosp_los ~0+ IcuAdmitYYYY:pf_cat + IcuAdmitYYYY + pf_cat + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),data = arf_hosp_los)arf_model_hosp_los_cat_tidy <-broom::tidy(arf_model_hosp_los_cat, conf.int = T)# Create predicted probabilities of hospital mortalityarf_hosp_los$pred_hosp_los <-predict(arf_model_hosp_los, type ="response")hosp_los_time <-ggplot(data=arf_hosp_los, aes(x=IcuAdmitYYYY, y=pred_hosp_los))+geom_smooth(data = arf_hosp_los, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf_hosp_los, aes(x=IcuAdmitYYYY, y=pred_hosp_los, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional Hospital LOS (Days)", title="Hospital Length of Stay over Time", color ="AHRF Category")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))hosp_los_arr <-arr(hosp_los_time)/100
Overall ICU length of stay is decreasing over time. Over the study period, there was a reduction in adjusted ICU LOS from 17.96259 (95% CI: 17.91665 - 18.00853) to 14.53623 (95% CI: 14.51806 - 14.5544). This was an absolute risk reduction (ARR) of 3.42636 (95% CI: 3.42586 - 3.42685). This change in outcome between the categories of AHRF over time were significant (95% CI did not cross 0 for interactions compared to None category of AHRF).
hosp_los_time
# Create DC Home Variablearf <- arf |>mutate(dc_home =if_else(hosp_dc_dest =="Home", 1,0))# Remove NAsarf_dc_home <- arf[complete.cases(arf$dc_home),]# Run MLMarf_model_dc_home <-glmer(dc_home ~0+ IcuAdmitYYYY + pf_cont + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),family =binomial(link ="logit"), data = arf_dc_home, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))arf_model_dc_home_tidy <-broom.mixed::tidy(arf_model_dc_home, conf.int = T, exponentiate=TRUE)# Run MLM for pg_cat interaction termarf_model_dc_home_cat <-glmer(dc_home ~0+ IcuAdmitYYYY:pf_cat + IcuAdmitYYYY + pf_cat + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),family =binomial(link ="logit"), data = arf_dc_home, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))arf_model_dc_home_cat_tidy <-broom.mixed::tidy(arf_model_dc_home_cat, conf.int = T, exponentiate=TRUE)# Create predicted probabilities of hospital mortalityarf_dc_home$pred_dc_home <-predict(arf_model_dc_home, type ="response")dc_home_time <-ggplot(data=arf_dc_home, aes(x=IcuAdmitYYYY, y=pred_dc_home))+geom_smooth(data = arf_dc_home, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf_dc_home, aes(x=IcuAdmitYYYY, y=pred_dc_home, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional Probability of Discharge Home", title="Discharge Home over Time", color ="AHRF Category")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))+coord_cartesian(ylim =c(0, 0.8))dc_home_arr <-arr(dc_home_time)
Overall discharge home is increasing over time. Over the study period, there was a reduction in adjusted hospital mortality from 66.744 (95% CI: 66.513 - 66.976) to 70.657 (95% CI: 70.569 - 70.746). This was an absolute risk reduction (ARR) of -3.913 (95% CI: -3.916 - -3.911). This change in outcome between the categories of AHRF over time were all not significant (all p>0.001 for interactions compared to None category of AHRF).
dc_home_time
# Create DC NH Variablearf <- arf |>mutate(dc_nh =if_else(hosp_dc_dest =="Nursing Home/Chronic Care/Palliative Care/Rehabilitation", 1,0))# Remove NAsarf_dc_nh <- arf[complete.cases(arf$dc_nh),]# Run MLMarf_model_dc_nh <-glmer(dc_nh ~0+ IcuAdmitYYYY + pf_cont + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),family =binomial(link ="logit"), data = arf_dc_nh, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))arf_model_dc_nh_tidy <-broom.mixed::tidy(arf_model_dc_nh, conf.int = T, exponentiate=TRUE)# Run MLM for pg_cat interaction termarf_model_dc_nh_cat <-glmer(dc_nh ~0+ IcuAdmitYYYY:pf_cat + IcuAdmitYYYY + pf_cat + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),family =binomial(link ="logit"), data = arf_dc_nh, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))arf_model_dc_nh_cat_tidy <-broom.mixed::tidy(arf_model_dc_nh_cat, conf.int = T, exponentiate=TRUE)# Create predicted probabilities of hospital mortalityarf_dc_nh$pred_dc_nh <-predict(arf_model_dc_nh, type ="response")dc_nh_time <-ggplot(data=arf_dc_nh, aes(x=IcuAdmitYYYY, y=pred_dc_nh))+geom_smooth(data = arf_dc_nh, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf_dc_nh, aes(x=IcuAdmitYYYY, y=pred_dc_nh, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional Probability of Discharge to Care", title="Discharge to Nursing Home/Palliative Care/Rehabilitation over Time", color ="AHRF Category")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))+coord_cartesian(ylim =c(0,0.2))dc_nh_arr <-arr(dc_nh_time)
Overall discharge to a facility is increasing over time. Over the study period, there was a reduction in adjusted hospital mortality from 8.973 (95% CI: 8.887 - 9.06) to 9.095 (95% CI: 9.062 - 9.128). This was an absolute risk reduction (ARR) of -0.122 (95% CI: -0.123 - -0.121). This change in outcome between the categories of AHRF over time were all not significant (all p>0.001 for interactions compared to None category of AHRF).
dc_nh_time
# Create DC Other Hospital Variablearf <- arf |>mutate(dc_other_hosp =if_else(hosp_dc_dest =="Other Acute Care Hospital", 1,0))# Remove NAsarf_dc_other_hosp <- arf[complete.cases(arf$dc_other_hosp),]# Run MLMarf_model_dc_other_hosp <-glmer(dc_other_hosp ~0+ IcuAdmitYYYY + pf_cont + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),family =binomial(link ="logit"), data = arf_dc_other_hosp, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))arf_model_dc_other_hosp_tidy <-broom.mixed::tidy(arf_model_dc_other_hosp, conf.int = T, exponentiate=TRUE)# Run MLM for pg_cat interaction termarf_model_dc_other_hosp_cat <-glmer(dc_other_hosp ~0+ IcuAdmitYYYY:pf_cat + IcuAdmitYYYY + pf_cat + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + (1|SiteID),family =binomial(link ="logit"), data = arf_dc_other_hosp, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))arf_model_dc_other_hosp_cat_tidy <-broom.mixed::tidy(arf_model_dc_other_hosp_cat, conf.int = T, exponentiate=TRUE)# Create predicted probabilities of hospital mortalityarf_dc_other_hosp$pred_dc_other_hosp <-predict(arf_model_dc_other_hosp, type ="response")dc_other_hosp_time <-ggplot(data=arf_dc_other_hosp, aes(x=IcuAdmitYYYY, y=pred_dc_other_hosp))+geom_smooth(data = arf_dc_other_hosp, aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(data = arf_dc_other_hosp, aes(x=IcuAdmitYYYY, y=pred_dc_other_hosp, color = pf_cat), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional Probability of Discharge to Other Hospital", title="Discharge to Other Acute Care Hospital over Time", color ="AHRF Category")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(4,1,1,1,1))))+scale_color_discrete(limits =c("Overall", "None","Mild", "Moderate", "Severe"))+coord_cartesian(ylim =c(0,0.2))dc_other_hosp_arr <-arr(dc_other_hosp_time)
Overall discharge to a different acute care hospital is decreasing over time. Over the study period, there was a reduction in adjusted hospital mortality from 11.47 (95% CI: 11.392 - 11.548) to 8.721 (95% CI: 8.691 - 8.751). This was an absolute risk reduction (ARR) of 2.749 (95% CI: 2.748 - 2.75). This change in outcome between the categories of AHRF over time were all not significant (all p<0.001 for interactions compared to None category of AHRF), with the exception of the severe category of AHRF.
dc_other_hosp_time
3.6.4 Secondary Outcomes: Time to Death (Adjusted Analysis) Truncated at 1 Year
arf |>filter(!(INV_IND=="NULL")) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = INV_IND))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="By Invasive Ventilation Status During ICU Admission", color ="Invasive Ventilation")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
arf |>filter(!(INV_DAYONE=="NULL")) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = INV_DAYONE))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="By Invasive Ventilation Status on Day 1 of ICU Admission", color ="Invasive Ventilation on Day 1")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
arf |>mutate(INV_AFTER_D1 =case_when( INV_DAYONE =="No"& INV_IND =="Yes"~"Yes", INV_DAYONE =="No"& INV_IND =="No"~"No")) |>drop_na(INV_AFTER_D1) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = INV_AFTER_D1))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="By Invasive Ventilation Status After Day 1 of ICU Admission", color ="Invasive Ventilation After Day 1")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
arf |>mutate(INV_TIME =case_when( INV_DAYONE =="No"& INV_IND =="Yes"~"After Day 1", INV_DAYONE =="Yes"~"Day 1")) |>drop_na(INV_TIME) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = INV_TIME))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="By Invasive Ventilation Status Timing", color ="Invasive Ventilation Timing")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
arf |>filter(!(NIV_IND=="NULL")) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = NIV_IND))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="By NIV Status During ICU Admission", color ="NIV")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
arf |>filter(!(ECMO_IND=="NULL")) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = ECMO_IND))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="By ECMO Status During ICU Admission", color ="ECMO")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
arf |>mutate(vent_sup =case_when( ECMO_IND =="Yes"~"ECMO", INV_IND =="Yes"~"IPPV", NIV_IND =="Yes"~"NIV",.default ="None")) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = vent_sup))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="By Ventilatory Support During ICU Admission", color ="Ventilatory Support")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_color_discrete(limits =c("None", "NIV","IPPV", "ECMO"))
# Figure 7- Mortality and Hospital Typeggplot(data=arf, aes(x=pf_cont, y=pred_freqmodel, color = HospitalClassification))+geom_smooth(data = arf, method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality",subtitle ="By Hospital Type", color ="Hospital Type")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
# Figure 8- Mortality and Genderggplot(data=arf, aes(x=pf_cont, y=pred_freqmodel, color = SEX))+geom_smooth(data = arf, method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality",subtitle ="By Gender", color ="Gender")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
arf |>filter(SEX==c("Male", "Female")) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = SEX))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality",subtitle ="By Sex", color ="Sex")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
Warning: There were 2 warnings in `filter()`.
The first warning was:
ℹ In argument: `SEX == c("Male", "Female")`.
Caused by warning in `==.default`:
! longer object length is not a multiple of shorter object length
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# Figure 9- Mortality and Agearf |>filter(!(age_cat =="NA")) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = age_cat))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality",subtitle ="By Age", color ="Age")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_color_discrete(limits =c("<44", "45-64","65-84", ">84"))
arf |>filter(!(ap3diag_min =="NA")) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = ap3diag_min))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality",subtitle ="By Admission Diagnosis", color ="Gender")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+scale_y_continuous(breaks =seq(0, 0.6, by =0.1))+scale_color_discrete(limits =c("Sepsis", "Medical","Trauma/Neurosurgery", "Post-Operative", "Cardiac Surgery"))
# Figure 10- Mortality and Treatment Limitationsarf |>mutate(TREAT_LMT =as.factor(TREAT_LMT)) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = TREAT_LMT))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="By Treatment Limiitation at Time of ICU Admission", color ="Treatment Limitation")+scale_color_manual(values =c("1"="lightblue", "2"="coral"),labels =c("1"="Full Active Management", "2"="Treatment Limitations"), na.translate =FALSE)+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
arf |>drop_na(EMG_RSP_ADM) |>ggplot(aes(x=pf_cont, y=pred_freqmodel, color = EMG_RSP_ADM))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality",subtitle ="By MET Call Admission", color ="MET Call Admission")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))
# Using IPPV at any timearf |>filter(!(INV_IND=="NULL")) |>ggplot(aes(x=IcuAdmitYYYY, y=pred_hosp_mort, color = INV_IND))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional Probability of Hospital Mortality", title="Mortality Trend of Patients Receiving IPPV Over Time", color ="Invasive Ventilation")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+coord_cartesian(ylim=c(0,0.2))
arf |>drop_na(NIV_IND) |>ggplot(aes(x=IcuAdmitYYYY, y=pred_hosp_mort, color = NIV_IND))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional Probability of Hospital Mortality", title="Mortality Trend of Patients Receiving NIV Over Time", color ="Non-Invasive Ventilation")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+coord_cartesian(ylim=c(0,0.40))
arf |>drop_na(ECMO_IND) |>ggplot(aes(x=IcuAdmitYYYY, y=pred_hosp_mort, color = ECMO_IND))+geom_smooth(method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional Probability of Hospital Mortality", title="Mortality Trend of Patients Receiving ECMO Over Time", color ="ECMO")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA)))+coord_cartesian(ylim=c(0,0.45))
arf |>mutate(vent_sup =case_when( ECMO_IND =="Yes"~"ECMO", INV_IND =="Yes"~"IPPV", NIV_IND =="Yes"~"NIV",.default ="None")) |>ggplot(aes(x=IcuAdmitYYYY, y=pred_hosp_mort))+geom_smooth(aes(color ="Overall"), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE, linetype =4)+geom_smooth(aes(x=IcuAdmitYYYY, y=pred_hosp_mort, color = vent_sup), method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x="Year", y="Conditional Probability of Hospital Mortality", title="Hospital Mortality over Time", color ="Ventilatory Support")+theme_bw() +theme(panel.grid.major.x =element_blank())+guides(color=guide_legend(override.aes=list(fill=NA, linetype =c(1,1,1,1,4))))+scale_color_discrete(limits =c("ECMO", "IPPV","NIV", "None", "Overall"))
ggplot(data=arf, aes(x=pf_cont, y=pred_freqmodel))+geom_smooth(data = arf, method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle="Using APACHE 3")+coord_cartesian(ylim =c(0,0.45))+theme_bw() +theme(panel.grid.major.x =element_blank())
# Run MLMarf_model_ap2 <-glmer(as.numeric(DIED_HOSP) ~0+ pf_cont + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + Apache2Score + (1|SiteID),family =binomial(link ="logit"), data = arf, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))# Create predicted probabilities of hospital mortalityarf$pred_ap2 <-predict(arf_model_ap2, type ="response")ggplot(data=arf, aes(x=pf_cont, y=pred_ap2))+geom_smooth(data = arf, method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="Using APACHE 2")+coord_cartesian(ylim =c(0,0.45))+theme_bw() +theme(panel.grid.major.x =element_blank())
# Run MLMarf_model_anzrod <-glmer(as.numeric(DIED_HOSP) ~0+ pf_cont + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + ANZRODRiskOfDeath + (1|SiteID),family =binomial(link ="logit"), data = arf, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))# Create predicted probabilities of hospital mortalityarf$pred_anzrod <-predict(arf_model_anzrod, type ="response")ggplot(data=arf, aes(x=pf_cont, y=pred_anzrod))+geom_smooth(data = arf, method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="Using ANZROD") +coord_cartesian(ylim =c(0,0.45))+theme_bw() +theme(panel.grid.major.x =element_blank())
# Run MLMarf_model_sofa <-glmer(as.numeric(DIED_HOSP) ~0+ pf_cont + SMOKINGINTENSITY + Apache3Score + CHR_RESP + CHR_CVS + SOFA_Score + (1|SiteID),family =binomial(link ="logit"), data = arf, nAGQ =0,glmerControl(optimizer ="bobyqa", optCtrl =list(maxfun =2e5)))# Create predicted probabilities of hospital mortalityarf$pred_sofa <-predict(arf_model_sofa, type ="response")ggplot(data=arf, aes(x=pf_cont, y=pred_sofa))+geom_smooth(data = arf, method ="glm", formula = y ~ splines::ns(x, 4), se=TRUE)+labs(x=expression(paste(PaO[2]:FiO[2], " Ratio")), y="Conditional Probability of Hospital Mortality", title="Association of Day 1 PF Ratio and Hospital Mortality", subtitle ="Using SOFA") +coord_cartesian(ylim =c(0,0.45))+theme_bw() +theme(panel.grid.major.x =element_blank())
3.10.2 Using Patients with > 12 Hours Invasive Ventilation
